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Eh ABSTRACT 

We investigate the global stability of a differentially rotating fluid shell threaded 
by vertical and azimuthal magnetic fields to linear, axisymmetric perturbations. This 
system, which models a thick accretion disk in the vicinity of its midplane, is susceptible 
to the Velikhov-Chandrasekhar (VC) instability in the absence of the azimuthal field. 
V^Q In most cases, the azimuthal field tends to stabilize the VC instability, although strong 

fields (Alfven speed of order the characteristic rotational speed in our incompressible 
model) are required for complete stabilization. Stability diagrams are constructed, indi- 
eating critical values of the two fields for instability. We find an additional strong field 
Qh instability that arises when the azimuthal Alfven speed exceeds the characteristic rota- 

tional speed. This instability, in the case of a freely bounded configuration, has certain 
^ similarities to the sausage instability for interpenetrating fields in plasma physics, and 

C3 may be important for very massive disks or filamentary molecular clouds. An applica- 

tion to the L1641 region in Orion A is briefly discussed. Finally, we find that the effect 
of a radially varying vertical field (without an azimuthal field) is mainly stabilizing. 
Subject headings: accretion, accretion disks - instabilities - ISM: magnetic 
fields - MHD 
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1. INTRODUCTION 

One of the more interesting recent developments in the theory of accretion disks was 
the discovery of virulent instabilities that develop only in the presence of magnetic fields. 
Balbus &l Hawley (1991) (hereafter BH) showed that a Keplerian disk in a state of pure 
rotation threaded by a weak axial field was subject to a local instability whose growth 
rate was on the order of the local rate of rotation. In a previous paper, we examined 
the global counterpart of this instability, the Velikhov-Chandrasekhar (VC) instability, 
showing that growth persisted at comparable rates (Curry, Pudritz, & Sutherland 1994, 
hereafter CPS). There remain serious questions, however, concerning how the instability 
is affected as models are augmented by additional physics. In particular, the influence 
of the more complicated magnetic field structures expected to exist in protostellar, CV, 
and AGN disks has yet to be carefully addressed. In this paper, we extend the model 
of CPS to include disks with radially varying vertical (B z ) and azimuthal (B^) fields. 

There are many reasons to expect an azimuthal field to be an important, sometimes 
dominant, magnetic field component in accretion disks. Strong differential rotation 
can generate B§ from a radial field component _B r , which can itself be created either 
by dynamo processes or accretion. The BH instability has been shown to generate 
strong B r and B<p from an initially weak B z (Hawley, Gammie, & Balbus 1995). If 
B z is inherited from the central object or the interstellar medium, disk torques can 
convert it directly to B^. Thus it is most likely that all three components of B are 
dynamically important for most types of disks. Of course, as has been made clear by 
all work following BH, one needs very little initial B z in order for that component to 
be "dynamically important." 

The observational evidence for B§ in protostellar disks is at the present time quite 
sparse, mainly due to uncertainties about the nature of the detected disks themselves. 
Recent mid-infrared spectropolarimetry of high-mass star-forming regions by Aitken et 
al. (1993) revealed a high correlation between objects with elongated molecular disk-like 
structures (numbering 10 in their sample) and magnetic fields oriented along the long 
axis of the disk (7 of these 10). The authors claim this as evidence for a predominantly 
azimuthal field structure in these regions. One should note, however, that the objects 
in question are 10 3 to 10 4 AU in extent, with masses ~ 10 3 Mq, and so are not likely 
to represent Keplerian accretion disks. As evidence for large-scale rotation is lacking, 
they may in fact be self-gravitating toroids or "pseudodisks," supported to some extent 
by the field itself (Galli & Shu 1993). 

Keplerian disks with magnetic field components in both the azimuthal and vertical 
directions have been actively studied as possible sources of centrifugally driven winds 
and outflows (e.g. Blandford Sz Payne 1982, Uchida Sz Shibata 1985, Pelletier Sz Pudritz 
1992). This suggests an additional motivation for the present study: to determine 
whether the various equilibria assumed in models of magnetically driven outflows are 
stable. A first step in this direction was taken recently by Lubow, Papaloizou & Pringle 



- 3 - 



(1994). 

In a different context, Galactic center molecular disk observations (Genzel 1989, 
Hilderbrand et al. 1990) indicate that B r ~ B§ ~ 1 mG, with a somewhat weaker 
B z . This is in contrast to the larger-scale field structure (i.e. the inner 70 pc of the 
Galaxy), which is almost purely vertical, i.e. perpendicular to the Galactic plane. 
Wardle & Konigl (1990) have modeled this region using a self-similar magnetized disk 
model, under the assumption that the inner field structure results from advection of the 
large-scale field by inflowing matter, with differential rotation subsequently leading to a 
strong azimuthal component. As observations of other galactic nuclei and AGN are still 
not able to resolve the inner disks, much less any associated magnetic field structure, it 
would be unwise to speculate further along these lines. However, since the inner regions 
of AGN are expected to possess "thick" rather than thin Keplerian disks, we use the 
same equilibrium sequence as introduced in CPS; namely, one in which radial pressure 
gradients oppose the central gravity for non-Keplerian rotation laws. The situation 
examined in the present paper is even more interesting, however, since radial magnetic 
gradients are also present. 

As a final possible application of the work presented here, we cite evidence that 
elongated filaments of gas in molecular clouds are associated with helical velocity and 
magnetic fields (Bally 1989). The latter are indicative of the simultaneous presence 
of and B z . The model employed in this paper, although formulated primarily for 
accretion disks, yields interesting results in the parameter range expected to hold in 
such regions. In particular, we find a new instability which sets in when the azimuthal 
Alfven speed is greater than the rotational speed. 

We defer to a later section a detailed description of previous work on the effect 
of B<p on the VC and BH instabilities, but it is of use to review here what is known 
generally about the stability of rotating configurations with azimuthal field. Since we 
do not attempt to account for the vertical structure of the disk in this study, the 
following discussion is restricted to purely radial distributions of angular velocity and 
magnetic field. The central question is this: given a rotation profile Q(r) and field 
distribution _B^(r), _B 2 (r), can one predict, even locally, whether a configuration is 
stable to infinitesimal perturbations in the fluid quantities ? What is needed is a 
necessary and sufficient criterion for stability, such as exists for purely vertical and 
purely azimuthal fields. These are, respectively; 

dfl 2 > 
dr ~ 

and 

1 d , o s9 r 2 d I Bs\ 2 
— (r 2 0) 2 -- r (-£ >0. (1.1) 



r 2 dr Aitp dr \ r 

The first criterion is due to Chandrasekhar (1960), and the second to Michael (1954). 
On the other hand, for the combined fields in the absence of rotation, the relevant 
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stability criterion is (Chandrasekhar 1961) 

£<r V < 0- 

In the most general case of combined fields with rotation, no similar criterion is 
known. Sufficient criteria are available, however; one of these is (Howard & Gupta 
1962, Dubrulle & Knobloch 1993, Kumar, Coleman, & Kley 1994): 

dr A-Kpr z dr 

Sufficient criteria can only be regarded as incomplete guides to the global stability of 
systems; the inherent limitations of criterion (1.2) will be made manifest later on in the 
paper. 

As to the actual distribution of and B z across the disk, there seem to be very 
few restrictions at this time 1 . By considering power-law distributions in these two field 
components, we hope to cover a range of plausibility. 

It is important to emphasize that the goal of this series of papers is not to replace 
the many local analyses that exist in the literature. Rather, a model such as we utilize 
below, while idealized and unrealistic in many respects, highlights intrinsically global 
behavior which will not be discovered in any local analysis. Examples of this found 
in the present work and in CPS are effects which involve coherent motions over large 
portions of the disk, and phenomena modified or even enhanced by the presence of a disk 
boundary, imperfectly modelled though it may be. Thus the present work complements, 
not replaces, existing local analyses. 

The format of the paper is as follows. The equilibrium state is described in §2, 
and the perturbations to this state in the following section. Quantitative results for 
the combined effect of azimuthal and constant vertical fields are presented in §4, and 
those for a radially varying vertical field in §5. In the final section, our results are 
compared with those of other investigators, and we make some additional comments on 
a new instability found in §4, before giving a final summary. Technical details of the 
calculations may be found in the four appendices. 

2. THE EQUILIBRIUM 

2.1. Basic Equations 

The equilibrium was described in detail in CPS, and has also been employed in 
stability analyses of thick, pressure-supported disks; see, e.g., Blaes & Glatzel (1986), 



1 Because an equilibrium B r immediately implies a time-dependent, growing B^, (BH) which destroys the time- 
invariance of the resulting equations, we ignore this field component in the present work. 
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Sekiya <fc Miyama (1988), and Jaroszyhski (1988). The model is a simplified form of 
the "thick torus" model for AGN (see, e.g., Paczyhsky &l Wiita 1980), supplemented by 
gradients of magnetic field pressure, but lacking vertical structure. It should therefore 
adequately describe a small region straddling the midplane of a real disk, with radial 
gas and magnetic pressure support taken fully into account. Thus the equilibrium is 
not that of a Keplerian disk, although this case is naturally included in the equilibrium 
sequence (CPS). 

Consider a cylindrical shell of homogeneous, incompressible, ideal MHD fluid, of 
infinite extent in the z-direction, rotating about the z-axis in the Newtonian point- 
mass potential \I/ = —GM/r. The purely radial dependence of the potential is justified 
if, at every radius r, the vertical scale height of the "disk" H <C r, so that there is 
little variation of \I/ with z. The stationary solution of the MHD equations depends 
only on the radial coordinate, r. To calculate explicit quantities of interest, we take 
the following power-law dependences for the angular velocity, azimuthal, and vertical 
magnetic fields, respectively: 

r \ — a / \ —6+1 / \ — c+1 



O(r) = O (^-J , B 4> (r) = B^[- ] j , B z (r) = B z0 (-J , (2.1) 

where Qq, B<po, B z o,a,b,c, and r$ are constants. As in CPS, we consider the effect of 
both rigid and free boundaries. In the latter case, B is supposed to permeate the regions 
both to the interior and exterior of the shell, as well as within the fluid. Using equations 
(2.1), the radial component of the equation of motion becomes (Appendix A) 



p - = ^ (-) - ™ - 1 (2 - w (-) + a - ^ (-) 

\r J r l r v \r J \r J 



1-26 / „ \ l-2c 



(2.2) 



where the prime symbol = d/dr, V^ z0 = B^ ^/Aitp are the azimuthal ((f)) and vertical 
(z) Alfven speeds at ro, and M is the mass of the central object (self-gravity is ignored). 

Inspection of equation (2.2) shows that the magnetic terms aid rotation and oppose 
the central gravity if b > 2 and/or c > 1, and vice- versa if b < 2 and/or c < 1. As in 
CPS, we consider configurations in which the gas pressure vanishes at the boundaries, 
and identify ro with the gas pressure maximum, where p' = 0. Equation (2.2) then 
gives 

9™. = ro Ql _ [(2 - 6 )y£ + (1 - c)V z l]/r . (2.3) 
r o 

This is merely a statement of radial magnetostatic equilibrium at the pressure maximum. 
In order for ro to be a maximum, we must have p"(ro) < 0. From equation (2.2), this 
requires 

(26 - 3)(6 - 2)Fj + (2c - 3)(c - 1)T^ + 2a - 3 > 0, (2.4) 

where an overline indicates that the Alfven speeds are now scaled with respect to ro^o, 
the circular speed at ro- Note that the above gives a > 3/2 when Vqq = and c = 1, as 
expected. Condition (2.4) should be satisfied for each equilibrium we examine. 
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Integrating equation (2.2) and eliminating GM via equation (2.3), one obtains the 
stationary pressure distribution 



pal 1 - r 

- = — + --1 + 7— 

p 2 r 2{a 



2(a-l) 



(2 - b)v; Q 



-(1 - c)V 2 



zO 



1 _ 1 - r- 2 ^- 1 ) 
r 1+ 2(6-1) 

] i _ r -2(c-i) 

- - 1 + ; ^~ 

r 2(c-l) 



(2.5) 



where we have chosen our units such that tq = Qq = 1, and where yu 2 is a constant 
equal to the ratio of thermal to kinetic energy at tq. We assume, as in CPS, that the 
gas and magnetic pressures remain finite as r<i — ► oo; this implies that a, 6, c > 1. 



2.2. Special Cases 

2.2.1. B z = constant 

Much of this paper is based on the particular case of a constant vertical field; i.e. 
c = 1. Then equation (2.3) and inequality (2.4) lead to the inequalities (in dimensionless 
units) 

1 - (2 - b)V^ > 0, 
(26-3)(6-2)T^ 2 + 2a-3>0, 

where we have dropped the overlines on the Alfven speeds for convenience. Considering 
all possible values of a and b leads to the conclusion that only certain values of V^q are 
permitted for a given (a, b). 

In the case of rigid boundaries, the inner and outer boundaries of the fluid are 
determined by the zeros of equation (2.5). For free boundaries, this is still true provided 
that B is continuous across the boundaries, and we shall assume that this is the case. 
Thus a given model is fixed by choosing r^/ri, a, b, and V^q. CPS found a monotonic 
increasing dependence of the VC instability growth rate on r^/ri, with a maximum at 
rijr\ J> 100 ; thus, we choose ^/n = 100 as a fiducial value for all calculations in this 
paper. The zeros of equation (2.5) can be positive, negative, or complex. The latter two 
(unacceptable) possibilities can occur even for (a, b, Vqq) obeying the above inequalities. 
We therefore conducted a three-parameter search for acceptable equilibria; the results 
are summarized in Figures la and b. 

In Fig. la, various critical values of the azimuthal Alfven speed are denoted by 
^1] V2, V3, . . .; each is a function of a and b. Although we calculated equilibria for 
all a, b in the range 1 < (a, b) < 3, Fig. lb shows only 3/2 < (a, b) < 2. We will 
restrict consideration for most of the paper to this range, since it reduces exactly to the 
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Fig. 1. — (a) Allowed regions and limiting azimuthal Alfven speeds in the (a, 6) plane, obtained 
from solution of the equilibrium equation (2.5) where p = 0, with rijr\ = 100. Equilibria for values 
of a and 6 lying in the shaded region and along dashed lines are not allowed for any Vqq. (b) 3D plot 
of allowed equilibria. Only the range 3/2 < (a, b) < 2 is shown. Permissible V^o f° r a given a, b lie 
between the upper and lower surfaces. The upper surface represents V§(a > b) and Ve(a < b). 

equilibrium of CPS when B<p = 0. The upper surface in Fig. lb represents V§(a > b) 
and Ve(a < b). Note that for a given a, 3/2 < a < 2, V5 > l^. 

There is another interesting property of the equilibrium relation (2.5): when 6=2 
and c = 1, we get the equilibrium of CPS. It may be checked that, for the power law 
fields we assume, this is the unique solution for which the current density, J = VxB/jdo, 
vanishes. Hence, the value of a is restricted to 3/2 < a < 2, just as in that study, and 
the location of the inner radius given by equation (2.9) of CPS. Since this special case 
allows us to examine the effect of the azimuthal field without the added complication 
of a current, we will assume 6=2, c = 1 (corresponding to ~ r _1 , B z = constant) 
when considering free boundaries in the sections to follow. 



2.2.2. B z = B z (r), B<p = 

In this case equations (2.3) and (2.4) give 

1 - (1 - c)V z 2 > 0, 
(2c - 3)(c - l)T/ 2 + 2a - 3 > 0. 

As above, these inequalities and equation (2.5) impose restrictions on allowed equilibria; 
these are summarized in Figures 2a and 2b. We now examine perturbations to the 
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above-described equilibria. 
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Fig. 2. — (a) Allowed regions and limiting vertical Alfven speeds in the (a, c) plane, (b) 3D plot of 
allowed equilibria. For the range of a shown, restrictions on V z o apply only for 1 < c < 3/2. The 
upper surface represents Vio- See text and Fig. 1 caption for details. 



3. THE PERTURBATIONS 

3.1. The Perturbation Equations 

We now consider the response of the above equilibrium state to small, axisymmetric, 
Eulerian perturbations of the form 

SX(r, z, t) = 6X(r)e i{kz+u}t \ (3.1) 

where X is any physical variable, and k and oj are the vertical wavenumber and frequency 
of the perturbation, respectively. Substituting the forms X + 6X along with equation 
(2.1) into the ideal MHD equations, linearizing, and eliminating all variables in favor of 
the radial velocity perturbation (see Appendix A for details), one obtains a second-order 
differential equation in Su r : 



where 



q(r) = k 2 r 



— [roj 2 (Su r )']' + q(r)8u r 



(3.2) 



r 2 



2V 



+ 



4k 2 (kVsVz 



Cj 2 



k 2 + 
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~2 _ 2 

UJ = UJ 



k 2 V 2 (r), 



and V^, ;Z (r) = B^ z (r) / \/4irp are the azimuthal and vertical Alfven speeds. The power- 
law form of the above (in dimensionless units) is 



q(r) = 2k 2 [bV$ 



-2b 



r — ar 



-2a 



+ (c-l)V z \ 



-2c 



4k 



+ —^-(kV^Vzo'. 



,1-b-c 



ujr 



-a\2 



UJ 



Cj 2 U 2 + 4 



(3.3) 



An alternative form of the perturbation equation useful for analytic purposes is 
obtained via the transformation 



whence equation (3.2) becomes 



i\) = {ruj 2 fl 2 6u r 



4>" = k 2 Q(r)4,, 



(3.4) 



with 
Q(r) 



UJ' 



{ar" 2a - bV 2 r~ 2b - (c - l)T/ 2 0? 



-2c 



UJ- 



-{kV^V^r 



l-b-c 



ujr 



+ 1 + 



fc 2 r 2 



1 



1 (r&) 



~,2y 2 



(ni 2 )" 



2 (ni 2 ) 2 



(3.5) 



When discussing free-boundary configurations, one must consider the form of the 
vacuum field perturbations in addition to those within the fluid. We will restrict our- 
selves to the current-free case, i.e. b = 2, c = 1, since then the perturbed magnetic field 
in the interior (r < r\, denoted by subscript i) and exterior (r > r-2, subscript o) regions 
is completely specified by a scalar potential such that 



<5Bj ;0 = B z Vxi,o-> 
Xi(zu) = cil (w), Xo(^) = c 2 K (w), 



(3.6) 



where \i,o = Xi,o( r ) ei ^ kz+wt \ w = l^l r i anQl an d Kq are modified Bessel functions of 
order zero. 



3.2. The Boundary Conditions 

We solve equation (3.2) subject to both rigid and free boundary conditions (BCs). 
The former are 

6u r (r\) = Su r (r2) = 0. 

Free BCs require the continuity of Lagrangian perturbations of the total normal stresses 
and magnetic flux across the boundaries. In the cylindrical geometry we are considering, 



- 10 - 



both B§ and B z are everywhere perpendicular to the surface normal n; thus B ■ n = 0, 
and provided that both fields are continuous across the boundaries, the appropriate BC 
is unchanged from the constant B z case; that is (CPS), 



(6u r )' + 



1 i ^ 2 ( i ; 2 T /2 Xi,o 

- + ~2 \ 9eff + k V z0 — 
X% o , 



r uj c 



Su r = 0. 



The subscript i applies at iq, subscript o at r 2 . 

For general power laws in B^r), B z (r), and in dimensionless units, the effective 
gravity is given by 



9eff 



= r 



-2 - ( 2 - f>)vy 



(1 - c)V z 2 r 



For vanishing current, this becomes identical to the g e fj of CPS; i.e. (? e // = r — 1/ 



From equation (3.6) one finds 

X 



1 I (wi) ^ Xo 



T~V\ 



r=T2 



1 K (w 2 ) 
\k\ Ki(vj 2 ) ' 



where w\ 2 = \k\r\ 2. 



4. RESULTS: CONSTANT VERTICAL FIELD 

The majority of our results have been obtained for the special case V z = constant. 



4.1. The Case of a = b 

When a = 6, the rotation frequency and its magnetic analog, V^,/r, have an 
identical scaling with radius. Equation (3.2) with (3.3) becomes 

^[r(6u r )']' + Q(r)6u r = 0, (4.1) 

where 

Q(r) = k 2 (Er- 2a -l)-±, (4.2) 

E = Jj [aCo 2 (V$ -l) + 2(n A V 4> -oj) 2 ] , (4.3) 

Oyl = kV z is the Alfven frequency, and where we have dropped the zero subscripts on 
and V z for convenience. The reader should note that these are constants throughout 
this section. 

Equation (4.1) with equation (4.2) is identical in form to the perturbation equation 
examined in CPS; for rigid BCs, the two problems are formally identical. The eigenvalue 
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E is a known function of a and k (see CPS and ff. equation (4.8)), but here its definition 
in terms of oj differs. The latter are solutions of the quartic polynomial obtained from 
equation (4.3): 

Eoj 4 - 2[EQ A + a(V$ - 1) + 2]oj 2 + m A V^oj + Vl\[EVl\ + 2a(V$ - 1) - 4V$] = 0. (4.4) 

The eigenvalue spectrum for E is infinite and every member is real and positive. For 
free BCs, the situation is complicated by the fact that oj appears in the BC itself. The 
problem is then no longer a standard Sturm-Liouville one and the introduction of E is 
not a particularly useful calculational device. In fact, when Y$ ^ 0, the resulting E is 
always complex. This leads to some interesting consequences which will be discussed 
presently. 2 For both sets of BCs, the resulting oj occur in complex conjugate pairs 
(Frieman & Rotenberg 1960). 

For simplicity, we begin by considering rigid BCs only. There are two special cases 
in which the roots of equation (4.4) have simple analytic forms. One is when = 1, 
which will be examined in the next section. The other is when a = 2. In that case, two 
of the roots are always stable, and the other two are 

w = ■ ( 4 - 5 ) 

The critical Alfven frequency for stability, where the imaginary part of equation (4.5) 
vanishes, is then 

1 ± (1 - TA 2 ) 1 / 2 

^A,crit = J-^j- 2 . (4.6) 

Two points are worth noting. First, has a stabilizing influence here. This could 
not be predicted from the sufficient criterion (1.2) given in the introduction, since b = 
2 => (rB^' = 0. Second, £lA,crit xs an explicit function of V^, and has two distinct non- 
zero solutions. That is, the stability criterion is altered in the presence of an azimuthal 
field, contrary to the claims of some recent investigators (see ff. §6.1). In the local limit, 
i.e. k — ► oo, V z — ► 0, we have E — ► r\ a = r\ (CPS); then equation (4.5) gives 



\[v^-(v 2 -2n A rl + n\rl)^] 



for the growing unstable mode. 

We solved equation (4.1) numerically, subject to both rigid and free boundary con- 
ditions, for a variety of (a, 6, V^) allowed by the equilibrium. As in CPS, we use the 
WKB approximation when V z <^ 0.3, since then the eigenfunctions 6u r are so sharply 
peaked that numerical solutions are difficult to obtain. 



2 When V4 = and for free BCs, as in CPS, it can be shown that the problem is still of Sturm-Liouville type, since 
to 2 is real and certain required conditions on the BC coefficients are satisfied; cf., Birkhoff & Rota (1989). 
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The principal results are as follows: 

(i) The VC instability persists for all V& < 1, 3/2 < a < 2, but with reduced growth 
rate. This conflicts with the naive prediction based on the sufficient criterion (1.2), 
since here (rB^)' > 0. The growth rate approaches zero as — ► 1. 

(ii) The presence of the azimuthal field also changes the stability criterion itself. Growth 
is damped at both short and long wavelengths. 

(iii) When > 1, a ^ 2, a new instability sets in, increasing in growth rate as V^. This 
large-field instability can be stabilized if V z is made sufficiently large. 

(iv) All of the unstable modes propagate; i.e. the real part of oj is ojr ~ kVj,V z . 

(v) The mode structure is unchanged from CPS; i.e. there exists a finite, ordered 
spectrum of unstable modes, whose growth rates are inversely proportional to a positive 
power of E. For the remainder of this paper, we will restrict consideration to the fastest- 
growing, or n = 0, mode. 




Fig. 3. — WKB growth rates as a function of Alfven frequency Qa = kV z for the fastest-growing 
(n = 0) mode, a range of V^, and two different vertical field values: (a) V z = 0.3; (b) V z = 0.05. 

In Figures 3a and b, we plot the dimensionless growth rate as a function of the 
Alfven frequency Qa = kV z . These curves show directly the effect of azimuthal field 
on the VC instability. We have chosen a = b = 2, but the curves are similar for other 
a = b. To display the effect for both strong and weak axial fields, Fig. 3a has V z = 0.3 
and Fig. 3b, V z = 0.05. Feature (i) is apparent in both figures; growth is clearly 
halted as 1^ -» 1. In the presence of _B^, growth rates are reduced due to vertical 
motions induced by magnetic pressure gradients (Blaes & Balbus 1994) (in the absence 
of B^, Su z ~ <5p; compare equation (3.1) of CPS and equation (A. 4), Appendix A). 
The instability couples to (stable) inertial modes, reducing its efficacy. The additional 
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stabilization provided by at shorter wavelengths (large 0^) is also apparent in both 
figures. The physical explanation for this is the same as in CPS; namely, that the 
restoring stress on a fluid element is more effective for distortions of larger curvature, 
i.e. at short wavelengths (also, see below). 

A new effect, the long-wavelength stabilization, is much more prominent in the weak 
axial field case (Fig. 3b). Even at ~ 0.7, one sees stabilization at long wavelengths 
(small Qa) for V z = 0.05. This behavior is entirely due to the presence of toroidal field 
lines, which provide an additional return force on a fluid element at long wavelengths. 
This can be seen by an explicit calculation of the perturbed magnetic tension, i.e. 

6(B ' V > B = ± (,kBJB - 2 SAi) , (4.7) 

where we have assumed without loss of generality that B§ ~ 1/r. As k — ► in the 
= case, the tension vanishes, indicating that instability persists up to the longest 
wavelengths. However, the second RHS term is independent of k, so that for nonzero 
B§ there exists an additional radial tension, which is always stabilizing. In addition, 
the effect is enhanced at small B z . It is this behavior that we observe at small in 
Fig. 3. 




Fig. 4. — Real (solid lines) and imaginary (dotted) parts of the eigenfrequency oj as a function of 
Alfven frequency for = 0.9 and (a) V z = 0.3; (b) V z = 0.05. Other parameters are the same as 
in Fig. 3. 

The real parts of all four roots for oj are shown in Figures 4a and b, for the same 
two values of V z and = 0.9. The corresponding imaginary parts are shown as dotted 
lines. The unstable modes (one growing, one damping) are created out of two real 
modes which merge for intermediate values of 0^. 
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Increasing to values in excess of 1 with a = b = 2 leads to no further instability. 
However, a new instability does occur for other values of a. The Keplerian case, for 
example, is shown in Figure 5. Each curve is labeled by its corresponding field values 
V<p, V z . As V z is held fixed at 0.3 and increased, the peak growth rate increases (solid 
curves). Were it not for the equilibrium constraint <^ 1.42 (see Fig. la), this growth 
would continue without bound as is increased. Now keeping fixed and increasing 
V z from 0.3 (dashed curves) leads to stabilization, until complete stability is achieved at 
V z ~ 0.81, implying (V z /V^) cr it — 0.57. We consider this large-field instability further 
in 84.3. 




Fig. 5. — Growth rates of the large-field instability as a function of Alfven frequency for Keplerian 
rotation and rigid BCs. Each curve is labelled by its corresponding V^, V z . Solid curves have 
V z = 0.3, while dashed curves have = 1.4. The chosen Alfven speeds are consistent with the 
equilibrium constraint < 1.42. 

Fig. 6. — Critical stability curves, Im oj = 0, in the (V^, V z ) plane for rigid BCs and (from right to 
left) a = b = 1.5, 1.7, 1.85, and 2. Growth rates increase from zero on both sides of each critical 
curve. The region at lower left is VC unstable for all a; the similar region at top left is large-field 
unstable for all a except a = 2. The dotted lines are the slopes of the LFI found analytically from 
equation (4.11). The large dots indicate upper limits on Y$ from equilibrium constraints; the curves 
are continued to larger for purposes of illustration. 



4.2. Critical Stability Curves 
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In CPS, it was shown that E behaves as 

E = ^^ + ■.■ + E (a), (4.8) 

where Eo(a) = T \ a an d E2(a) = lim^o^ 2 E . Since the longest wavelength perturbations 
are always unstable in that case, one could then calculate the critical field strength for 
stability, V Z}Cr it, by taking the limit of the dispersion relation as oj — ► 0, k — ► 0. When 
V<f, ^ 0, the values of (V^, V z ) for which marginal stability holds constitute curves in the 
(V^, V z ) plane. This section will be concerned with the construction of such curves. 

When Vtj, ^ 0, there is an added complication. Fig. 3 shows that the most persistent 
unstable mode is not always that with k — ► 0. Rather, the last unstable mode which 
persists as V^ — ► 1 has intermediate k\ the precise value is a function of V z . We note 
here that in the local limit, k — ► oo, V z — ► 0, the growth rate curves are perfectly 
symmetrical about Q A = 4, which is the value for peak growth when a = 2 and V$ = 0. 

When V§ = 1, equation (4.4) has the four roots 

n A (twice), ±2E~ 1 / 2 - Q A , (4.9) 

all of which are real, provided that E is real. This result is independent of both V z and 
k. Thus the line = 1 must lie in an absolutely stable region in the (V^, V z ) plane. 
Further, taking = 1 ± e with e small and positive, and expanding oj in orders of e, 
one finds from the first-order correction that as long as a < 2, instability occurs. Thus 
the line V§ = 1 constitutes an absolutely stable region in the (V^, V z ) plane. We need 
not take any special care when considering V z — ► 0. 

For larger values of V z , say V z J> 0.3, the limit k — > does give a reliable estimate 
of the critical curve (Fig. 3a). Now taking 

oj = koj\ + k^OJi + fc 3 W3 + . . . 

along with equation (4.8), equation (4.4) becomes, to first order, 

E 2 oj\ -2[E 2 V Z 2 + a(Vl - 1) + 2]oj\ + SV^^ + V z [EV Z + 2a(V$ - 1) - 4V$] = 0. (4.10) 

Solving equation (4.10) for the loci of oj\ =0 in the (V^,, V z ) plane gives the critical 
stability curves we seek. These have been plotted in Figure 6. The = results, 
which were derived in CPS, are obtained where curves intersect the V z axis. One sees 
that as Vfj, is increased from zero, smaller values of V z are needed for stabilization, until 
at Vfj, = 1 complete stabilization occurs for all a. For a = 2, all V^ > 1 are completely 
stable; this is represented by the heavy line along the a = 2 curve and continuing along 
the V§ axis from V§ = 1 to infinity. For a < 2, the plane above V$ = 1 is divided into an 
unstable part (adjoining the axis) and a stable part (to the right of a given critical 
curve). The unstable region at > 1 extends to infinity and shrinks to zero size as 
a — ► 2. Actually, the size of the unstable region for a ^ 2 depends on the particular 
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value of a (and when a ^ 6, on b as well). This is due to the equilibrium constraints 
placed on Y$ by Fig. 2b. The largest allowed for each a has been indicated in Fig. 6 
by a large dot on the appropriate critical curve. We extend the curves to higher values 
of Vtj, merely to display their asymptotic behavior (see below); such large field values 
will not be attainable in reality. 

4.3. The Large-Field Instability 

The almost linear behavior of the curves in Fig. 6 at large V^, V z is intriguing. In 
this limit, and again taking k — ► 0, equation (4.4) gives 



.2 



Wi 



n 2 A [E 2 n 2 A + 2( a - 2)n 2 



1 2P(E 2 n 2 A + aSll) ' 

where 0^, = kV<p. This implies 

V z /V^ > (V z /V^) cnt = ^2(2 - a)/E 2 (4.11) 

for stability. Values of E 2 and iV z /V^) CT it for 3/2 < a < 2 may be found in Table 1. 
The equality in (4.11) gives the asymptotic (i.e. large V<p,V z ) behavior of the critical 
stability curves, as shown by the dotted lines in Fig. 6. 

The nature of this large-field instability (LFI) is easily understood upon comparison 
with the equivalent nonrotating system. The equilibrium pressure distributions are 
compared in Appendix B, where it is shown that when V§ ^> foOo, the system reduces 
to its nonrotating equivalent. For the latter, Chandrasekhar (1961) derived the necessary 
and sufficient stability criterion 



a 


E 2 (a) 




1.5 


2.55 


0.63 


1.6 


2.78 


0.54 


1.7 


3.01 


0.45 


1.8 


3.26 


0.35 


1.9 


3.50 


0.24 


2.0 


3.75 


0. 



Table 1: Ratio of critical Alfven speeds, (V z /V^) cr it : as a function of shear parameter a for the LFI. 
See text for the definition of E 2 (a). 
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where I\ is a positive-definite integral function of r. 3 Differentiating both sides of this 
inequality, and assuming that k ^> d/dr (this is equivalent to considering the longest- 
wavelength radial perturbations, which should be the most unstable), we obtain its local 
version; 

k2B z > ^T-(^) 2 = — J z> ( 4 - 13 ) 

where J z is the axial current. The LHS of (4.13) represents the restoring force exerted 
on a radially displaced fluid element by the perturbed vertical field, while the RHS 
is the excess Lorentz force on that element due to perturbations of B§. The latter 
is the exact analogue of the destabilizing centrifugal force in the BE instability. Since 
J z = (2 — b)B^,/r, configurations with ^> tqVIq and b > 2 are stable to the LFI. 
In essence, the LFI is the result of an imbalance between radial gravity and a radially 
stratified, buoyant magnetic field (see also Appendix B). 

4.4. Free Boundaries 

The only case to be considered here is a = b = 2, since we restrict consideration 
to the current-free situation. The critical stability curve is shown in Figure 7. The 
most significant difference is the disappearance of the absolutely stable line at V<p=l. A 
glance back at the roots (4.9) of the polynomial (4.4) shows how this happens. When 
V<f, ^ 0, E is no longer real, and one of these roots becomes growing unstable. The 
actual behavior is as follows. Consider a line of constant V z , such that < V z < 1. 
The peak growth rate for a given decreases from a maximum at = 0, to some 
minimum in the vicinity of ~ 1, and then increases again without bound as is 
made larger. Note also how much more extended is the unstable region in the free case 
versus the rigid one. 

Global effects must clearly be at work here, since Va, > 1 is unstable only in the free- 
boundary case. As rotation is not likely to be important in this region, it is instructive to 
consider the equivalent nonrotating problem. A situation similar, although not identical, 
to the latter is that of the plasma "pinch" (e.g., Chandrasekhar 1961, Ch. XII, §115). 
This consists of a filled cylindrical column of plasma, threaded by a uniform B Zl and 
surrounded by a vacuum region containing the same B z together with an azimuthal 
field oc r _1 . The entire arrangement is usually encircled by a concentric conducting 
wall, but we are free to place this at infinity and so ignore it for our present purpose. 



3 Compare criterion (4.12) with that for the VC instability; i.e. 
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V(r„ ) V 2 /(r n ) 



Fig. 7. — Same as Fig. 6, but for free BCs. From lower right to top left, curves are for a = 
1.5, 1.7, 1.85, and 2. The region to the left of each curve is unstable; that to the right, stable. 

Fig. 8. — Same as Fig. 6, but for a ^ b. Each curve is labelled by its corresponding a, b. For a > b, 
unstable regions lie to the left of each curve, stable regions to the right. For each a < 6, there are 
two branches of the critical stability curve. One, at Y$ < 1, bounds the VC unstable region from 
above; the other, at > 1, bounds the large-field unstable region from below. The large dots 
indicate upper limits on from equilibrium constraints. 

For the extended configurations we consider (r^/ri = 100), the (nonrotating) situation 
is nearly identical except for the fact that in our problem B z and B<p interpenetrate 
everywhere, not just in the vacuum region. However, such interpenetrating fields have 
been considered by Tayler (1957), with the finding that such arrangements are more 
unstable. 

When all fields are continuous across the plasma/vacuum boundary, the fluid is 
susceptible to the well-known (m = 0) sausage instability, which can be stabilized if 
and only if V z > V?/2 =4* (V z /VA))crit <^ 0.707. It is of interest to compare this 
figure with the inverse of the slope of the critical curve for a = b = 2 in Fig. 7, 
which is (V z /V^) cr it ~ 1.5. The latter situation is more unstable, we posit, due to the 
interpenetration of B<p and B z in the fluid region. Since the exterior B<p is the cause 
of the sausage instability in the first place, it is not hard to imagine that its presence 
inside the fluid will inhibit the stabilizing effect of B z . 



4.5. The General Case: a ^ b 
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4-5.1. Rigid Boundaries 



When a ^ b, the reduction of the full eigenvalue problem, equations (3.2) and (3.3), 
to a single characteristic polynomial is no longer possible. Before proceeding to a nu- 
merical solution, however, it is of use to present such analytic formulae as are available. 
There are two approaches which have had some success in this regard, and which lead 
to identical results. One is the local analysis of Dubrulle & Knobloch (1993), which 
ignores radial variations in equilibrium quantities compared with those of perturbed 
ones (i.e. r(6X)'/6X 1). The other is the slender annulus approximation adopted by 
Kumar, Coleman, <fc Kley (1994) which we follow here to preserve the global character 
of the analysis (Appendix C). 4 In the limit V z — ► 0, both give the following condition 
for stability, 

[2-a-(2-b)V*](a-bV>)<0. (4.14) 
For example, if < b < 2 and b < a, stability holds if 

— b < v i< v (4 - 15) 

whereas for b in the same range and a < b, stability holds if 

a n 2 — a , 
j<V}< — b . (4.16) 

It is easy to see that both of these inequalities bracket = 1. 

As regards the (V^,, V z ) critical stability plane, equations (4.15) and (4.16) imply 
the existence of a stable region along the axis bracketing = 1. How this limiting 
behavior is related to the critical curves for general V^, V z , and rijr\ will now be 
investigated. 

For configurations with rigid boundaries, all (a, 6, V^) consistent with Fig. lb may 
be considered. Qualitatively, there are some significant differences from the a = b case. 
These differences may be classified according as a > b or a < b. Several representative 
critical stability curves are shown in Figure 8. Beginning at the far right-hand side of 
the diagram, we have a stable region at large V z . When a > b, the curves achieve a 
minimum value of V z for some J> 1, and then display the linear asymptotic behavior 
found in the previous section. For a > b, there exists no stable region (not even = 1 ) 
at small V z . This result contradicts the local prediction (4.15) of a stable region as 
V z — ► 0. As a is reduced to values nearer to 6, e.g. a = 1.7, b = 1.6, the "knee" of the 
curve bends inward to smaller values of V z \ it is easy to imagine what happens in the 
limit as a — ► b from above; the knee of the a > b curve deforms into the line V<p = 1, 
which extends all the way to V z = as in Fig. 6. 



4 This approach is actually superficially global, in that although radial BCs are applied, the authors assume that 
the boundary separation is proportional to \/k, k ^> 1. 
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If a is decreased further such that a < b, the situation is less clear. We have been 
able to confirm numerically the persistence of two distinct unstable regions, one at 

J> 1 (LF unstable) and one at <^ 1 (VC unstable), down to values of V z ~ 0.2. 
Between the two stability curves lies an absolutely stable region, bracketing = 1. At 
smaller 14, mode crossing becomes a significant hindrance to the numerical algorithm, 
and precise determination of the critical curves is difficult. For a = 1.7, b = 1.75, we 
were able to follow the n = mode down to V z ~ 0.2 (solid curves in Fig. 8); beyond 
this, we join the numerical curves onto the values given by the local relation (4.16) at 
V z = (dashed curves). 

It should be mentioned that this region of parameter space, i.e., 

V z -+ 0, « 1, a < b, 

is highly restricted by the equilibrium constraints. A glance at Fig. la reveals that we 
must have a > 3/2. Since the LFI requires b < 2, we therefore have 3/2 < a < 2, a < b 
as our region of interest. Widely separated values of a and b in this range have limiting 
Alfven speeds well below unity; e.g. when a = 1.55, b= 1.95, Ve = 0.46. Hence, the LFI 
is not a concern. Less separated values of the two parameters allow larger equilibrium 
fields; e.g., a = 1.85, b = 1.95 =4* Vq = 2.93. But it is likely that for such a, b the 
critical stability curves are qualitatively similar to the a = 1.7, b = 1.75 case shown in 
Fig. 8. To confirm this, we developed an approximation whose validity depends on the 
smallness of the parameter a/b, but imposes no restrictions whatsoever on the global 
geometry. 5 The critical stability curves found by this method always contain a stable 
region bracketing V§ = 1. 

To explain the existence of a stable region at small V z , it is instructive to look at 
the dependence of the perturbations on a and b. The VC instability arises from an 
imbalance of the destabilizing stress B z 8B§j '4ir and the stabilizing stress B z 6B r / '4ir . 
When an azimuthal field is present, the ratio of these as found from equations (A. 7) 
and (A. 8) is 

8 1± = \ 2 kB z Q - ^(bu 2 + 2Q2 )1 . (4.17) 

6B r kB z u 2 [ z ruj K AJ \ V ' 

The first term in the square brackets behaves as r~ a , the second term as r~ b . Consider 
unstable modes only, so that oj ~ (this is still true when <^ 1). The relative 
magnitude of the two terms then depends on: (a) the relative magnitude of a and b, (b) 
the relative magnitude of and B^/r, and (c) whether r < 1 orr > 1 (i.e. inside or 
outside the pressure maximum, respectively). Assume that J> B^/r; i.e. that we are 
in the VC regime. Recall from CPS how strongly peaked were the radial eigenfunctions 



6 Specifically, we define a new variable, x = 1 — W 6-1 , a < b, and expand the perturbation equation (3.4) in 
powers of x. Finding a series solution and subjecting it to rigid BCs, one obtains a fourth-order dispersion relation 
similar to equation (4.10), which can be solved numerically for to. See Curry (1995) for details. 
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of the unstable modes interior to the pressure maximum; this suggests that the region 
r < 1 is far more important than r > 1 for the linear stage of instability. We therefore 
restrict consideration to that region. Now, when a > 6, the first term in the above 
dominates the second, and SB^/SBr retains the same sign as it had in the absence of 
B^, where its effects were always destabilizing. Thus while one would expect a reduction 
in the growth rate near ~ 1, it should not completely vanish. 

On the other hand, when a < b, the second term in equation (4.17) can be compa- 
rable to the first even when B^/r <^ 0, and so a change of sign in SB^/SBj. occurs at 
some Va <^ 1, signifying stabilization. Such stabilization cannot be maintained at higher 
values of V^, however, once the LFI begins to set in. This gives the upper boundary of 
the stable region. Physically, B^ overwhelms in the inner disk when b > a, leading to 
momentary stabilization until the field becomes so strong that rotation is no longer a 
viable means of support. At this point, the LFI takes over. Because the local analysis 
gives no information about the radial dependence of the eigenfunctions, Dubrulle & 
Knobloch were not able to detect this interesting dependence of the stability properties 
on the relative magnitudes of a and b. 

4-5.2. Inapplicability of the Local Approximation 

There are two particular cases in which the local criterion predicts qualitatively dif- 
ferent behavior than that examined above. When either a = 2 or b = 2, criterion (4.14) 
yields the following results: 

(i) a = 2, b > 2; > ^/2/b, 

(ii) a = 2, b < 2; V<j, < y/5Jb, 

(iii) b = 2, a < 2; V<p > y/a/2, 

(iv) b = 2, a > 2; V<p < ^fi~/2. 

In all these cases, there exists only a single critical curve. Since (ii) and (iv) have 
a > 6, we expect the local prediction to be unreliable by extension of the results of the 
previous section; thus we do not expect a critical stability curve to extend all the way 
to V z = 0. In cases (i) and (iii), however, there is no a priori reason to doubt the local 
results. 

As test cases, consider the physically interesting power law indices 

(i) a = 2, b = 3; V<p > 0.82, 
(iii) 6 = 2, a = 1.5; ^Vj > 0.87. 

The first case is that of constant angular momentum, with a rapidly decreasing az- 
imuthal field. The second is a zero-current, Keplerian configuration. By the results of 



- 22 - 



§4.3., both systems should be stable to the LFI. In addition, both the Michael (equation 
(1.1)) and the Howard & Gupta (equation (1-2)) criteria are satisfied. The equilibrium 
constraints place no restrictions on the value of Ya for these (a, 6). The critical stability 
curves, calculated numerically, are shown in Figure 9. Again, it is difficult to extend 
the curves much past V z <^ 0.2, but in case (i) we have been able (quite remarkably) 
to follow the curve down to V z = 0.05. As in the rest of the paper, the results are for 
n = 0, which we have always found to be the fastest growing radial mode. 




V(r„ n„) 



Fig. 9. — Critical stability curves for two special cases examined in §4.5.2. The dashed and dot- 
ted lines indicate the local predictions; they intersect the V<p axis at Y$ = 0.87 and Y$ = 0.82, 
respectively 

Fig. 10. — Selected eigenfunctions at peak growth for a = 2, 6 = 3, and increasing V<p from top left 
to lower right. The solid line indicates the real part of Su r ; the dashed line, the imaginary part. 
Each eigenfunction is normalized to its peak value. 

The results are surprising in that they bear no resemblance to the local predictions 
(the dashed and dotted lines in Fig. 9) as V z — ► 0. The a = 2, 6=3 curve, e.g., shows 
that the VC unstable region is five times as large at V z = 0.1 than the local prediction, 
and the curve even appears to be diverging as V z — ► 0, instead of approaching a constant 
value. One reason why the local approximation fails here can be found via inspection of 
the relevant eigenfunctions, a few of which are plotted in Figure 10. When Y$ ^ 0, Su r 
has both real (solid line) and complex (dashed line) components. As V& is increased from 
zero, one sees a gradual spreading of the eigenfunction from the inside regions outward. 
In the region near the critical curves as V z — ► 0, Su r is much more extended than in any 
other case examined thus far. The peak of the eigenfunction at maximum growth is no 
longer confined to the small region between r\ and ro; e.g. when V z = 0.3, = 1.5 
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(Fig. lOd), it lies at v/vq m 3, and 6u r has a nonnegligible amplitude over the entire 
shell. This feature alone is enough to show that the local and thin shell analyses are 
inadequate to capture the true behavior of the system in this parameter regime. It also 
confirms one of the main findings of CPS; namely, that the local and 'critical' limits are 
antipodal: the latter can only be reached via a global analysis. 

4-5.3. Free Boundaries 

For 6=2, we plot a variety of a values in Fig. 7. Again, in contrast to the rigid BC 
case, there exists no stable region around V& = 1; this is easily understood in light of 
the discussion given in §4.4. The unstable regions are larger for a ^ 2 than for a = 2; 
this is due to the fact that two instabilities, the current-driven LFI, and the sausage 
instability, act simultaneously. The asymptotic critical values for these curves range 
from (Vz/V^crit = 1.5 for a = 2 to (V z /V^) cr it = 2.3 for a = 1.5. The l^-axis intercepts 
of the curves match the values found in CPS. 

4.6. The Effect of Simulated Vertical Boundaries 

In CPS, we calculated the critical V z for several fixed, nonzero values of k, corre- 
sponding to vertical wavelengths, \ cr it, between 100 and 0.1 in units of the inner radius 
r\. The intent was to gauge the probable effect of vertical disk boundaries on V Z criU 
under the hypothesis that the longest unstable wavelength could not exceed the disk 
thickness. The interesting result was that for \ C rit = 0.1, a reasonable value for a thin 
Keplerian disk, V ZjCr it — 0.04 r=i V z j{, where V z j{ = "v/6 c s /ir is the local Keplerian crit- 
ical field estimate (BH). Thus, the super-rotational Alfven speed required for stability 
in the infinite incompressible cylindrical shell model translates to a super-thermal V z in 
a thin, isothermal disk. 

In the presence of an azimuthal field, we have found that for small V z , values of 
V(j> ~ r 0^0 are required for critical stability. This therefore begs the same question as 
asked in CPS: does the same result hold for thin disks, or does critical stability again 
require ~ c s ? 

Following the same calculational procedure as in CPS, we calculated critical stability 
curves for a = b = 2, rigid boundaries, and a range of \ cr it (Figure 11). Mode confusion 
prevents us from going to \ C r%t < 0.2, but the trend is clear. The curves do not all 
approach = 1 as V z — ► 0, since they are for fixed k\ the small Qa stabilization 
discussed in §4.1. takes over when V z becomes small. This can be seen explicitly by 
deriving the following "local" critical stability relation. In the local limit, E — ► rf, so 
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E l l 2 a (0.5) 2 = 0.25 for a = 2, and equation (4.6) gives (in proper units) 



l crit 
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Assuming the azimuthal field is subthermal so that it does not significantly alter the 
overall structure of the disk, the critical stability requirement \ cr it ~ 2iJ = 2v / 2c s /Q(ri) 
then yields 
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(4.18) 



For \ cr it = 0.1 iq, equation (4.18) gives the long-dashed curve shown in Figure 11. 
Although equation (4.18) concurs with the sequence of curves shown and highlights 
their key qualitative features, it cannot be rigorously correct for two reasons: first, one 
cannot actually have a "thin" disk with a = 2; and second, the derivation is inconsistent 
for V^/r Qo ~ 1 > V^/cs. 




Fig. 11. — Critical stability curves for selected perturbation wavelengths A cr ^, for a = b = 2 
and rigid BCs. The long-dashed curve is the "local" critical curve given by equation (4.18) with 
\ cr it = 0.1 ri. 

Fig. 12. — Numerical growth rates as a function of Alfven frequency for = and V z = 
V z or 1 ~ c (n = mode). Curves are labelled by their corresponding c values. 

Irregardless of the applicability of equation (4.18), the numerical curves in Figure 11 
unambiguously show that although V ZjCr it decreases with decreasing \ C rit ( or decreasing 
scale height H), the same is not true of V^^rit- Even in the thin disk limit, one still 
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requires V^^rit ~ r 0^0 f° r complete stabilization; i.e. for all wavelengths and at any V z . 
This result can be understood by recalling the physical cause of the LFI: it can only 
occur when rotation is relatively unimportant in comparison with the azimuthal field, 
a requirement that does not change when the effective scale height is reduced. 

In a real, compressible, vertically stratified accretion disk, Parker (vertical magnetic 
buoyancy) instability is known to act when J> c s <C ^oOo- Thus, the above result 
could have at least two important consequences for such a disk. First, it argues per- 
suasively against the possibility of the LFI ever occurring, since for J> c s , Parker 
instability would already have caused a rearrangement of the magnetic equilibrium. 
Second, and more importantly, the above result suggests that the VC instability is un- 
likely to be stabilized by an azimuthal field of any power-law index or strength <^ c s . 
We will discuss other possible environments for the LFI in §6.2. 

5. NONCONSTANT VERTICAL FIELD 

Should an accretion disk be threaded by a vertical magnetic field, the latter is more 
likely to vary with radius than be uniform. Although we do not explicitly model the 
accretion flow in this study, its overall effect is to drag field lines radially inward (by 
flux-freezing), leading to a higher B z flux in the inner regions. In this section we consider 
the effect of a radially varying vertical magnetic field on the VC instability, and neglect 
the azimuthal field. Although for completeness it would be desirable to consider the 
most general situation of nonconstant vertical and azimuthal fields, we defer that to 
a future work. An additional complication arises in that case, since resonances can 
occur where the real part of uj 2 — k 2 V 2 (r) = 0. This is not a concern for the unstable 
modes considered in this section, since they always have uj 2 < 0; a proof of this is given 
in Appendix D. We consider only rigid BCs, since the zero-current restriction on our 
freely-bounded equilibria requires c = 1. 

Even with the restriction to rigid BCs and the knowledge that uj 2 is real, analytic 
progress is difficult, since the r-dependence of uj 2 means that the perturbation equation 
(3.2) is not of standard Sturm-Lioville type. Regrettably, the global WKB approach 
used in CPS does not give satisfactory results in this case, for the following reason. 
Choosing 1/k as a small parameter, the last RHS term of equation (3.5) cannot be 
neglected, due to the presence of w-dependent terms. However, use of the thin shell 
approximation of Appendix C gives the result 

2 Oy2(fc 2 + c) + r]/2] + 2(2 - a)k 2 
" ~ 2(fc 2 + 3/4) 

± i^Ao[k 2 (v + 1) + Viv/* + 2c " 3/4) + 3(1 - 2c) + 4c 2 ] 

+ 2k 2 n 2 A0 [8k 2 + 3a + (2 - a)(4c + r,)} + 4k 4 (2 - a) 2 } 1 ' 2 /2(k 2 + 3/4), 

where -q = 3 — 8c + 4c 2 and Qao = kV z (ro). This solution can be regarded as quan- 
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titatively valid only in a small neighbourhood of the pressure maximum. However, it 
exhibits roughly the same qualitative behavior as the exact numerical solutions dis- 
cussed below. In addition, taking k 3> 1 leads to the local dispersion relation of CPS 
and BH. 

Exact numerical growth rates as a function of O^o f° r various values of c > 1 and 
the fiducial values a = 2, V z o = 0.3, 1*2/1*1 = 100 are plotted in Figure 12. The different 
curves are labelled by their corresponding c values. For c > 1, the growth rate is always 
reduced from its constant V z value. The critical Alfven frequency for stability, ViA0,crit, 
decreases with increasing c, until at some critical value, c 2.5 — 3, it begins to increase 
again. The peak growth rate, however, continues to decrease. We have difficulty finding 
\uj\ for c J> 3.5 and large Q^Oi possibly due to the simultaneous presence of several 
unstable modes with the same growth rate. 

The particular laws c = 9/4 and c = 5/2 correspond to the flux distributions for two 
popular centrifugally-driven wind models; Blandford & Payne (1982) and Pelletier & 
Pudritz (1992), respectively. As far as the stability of these distributions is concerned, 
there is no great distinction between either; both are VC unstable. One should note, 
however, that both models require B<p ^ 0; in the former B<p ~ B z , while in the latter, 

~ r~ l . Thus while the results of the present paper suggest that V<p ^ i*o^o will 
further stabilize, a calculation explicitly incorporating B^ is still necessary. 

The run of peak growth rate with c is shown in Figure 13, for V z = 0.3, 1*2/1*1 = 100, 
and different values of a. The a = 1.5 curve is incomplete because 1 < c < 3/2 is 
forbidden by the equilibrium (Fig. 2b). The large dot on the vertical axis shows the 
constant V z value (see Fig. 7a of CPS). 

The physical reason for the stabilization observed here is the same as for an azimuthal 
field in the presence of a constant V z , except that now the additional vertical motions are 
induced by the gradient of the vertical field (see equation (A. 4)). As for the effect on the 
stability criterion, we advance the following argument. In the inner region of the disk, 
V z (r < 1*0) > T4(i*o)- Thus, the local instability at r < r$ will be attenuated compared 
to the constant V z = V z (tq) situation. By the same argument, the local growth rate 
should be enhanced outside the pressure maximum. However, the unstable eigenmodes 
are strongly peaked inside r = 1*0 1 when c « 1; this region is more important for the 
action of the VC instability. Thus for c J> 1, the attenuation effect dominates, and 
the critical wavenumber for stability, k cr n = CIao crit/VzOi is reduced from its value in 
CPS. For c significantly greater than 1, an interesting phenomenon occurs (Figure 14). 
The peak of the eigenfunction Su r gradually moves from inside the pressure maximum 
(for c fa 1) to r/i*o ~ 1 (for c = 7/2), and presumably beyond for more extreme field 
gradients. Thus it is likely that for larger c, the above argument no longer holds. That 
is, the enhancement effect of the VC instability at r > r$ does contribute, leading to a 
reversal in the trend of CIao crit- 

We have searched for other unstable modes, e.g. at V z ^> 1, with no success. 
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c r/r 

Fig. 13. — Maximum growth rates as a function of vertical field index c for different rotation indices 
a. The lower curve is incomplete since equilibria with a = 1.5 and 1 < c < 1.5 are not allowed (Fig. 
2a). The maximum growth rate at c = 1, found in CPS, is shown as a large dot slightly offset (for 
clarity) from the vertical axis. 

Fig. 14. — Radial eigenfunctions (n = mode) at peak growth rate for various c. The overall 
normalizations have been adjusted to unity for purposes of comparison. 

Interchange modes, which might be expected to act at large field strengths, do not occur 
here because we consider only axisymmetric perturbations (see, e.g., Kaisig, Tajima, & 
Lovelace 1992, Lubow <fc Spruit 1995). 

6. DISCUSSION 
6.1. Comparison with Previous Results 

Here we compare our results for the effect of the azimuthal field on the VC instability 
with those of four recent papers, finding some significant discrepancies. 

Dubrulle Sz Knobloch (1993) (DK), via a WKB method, found that the imaginary 
part of the eigenfrequency, ojj ~ Q A /(1+ const, x V%) in the limit V z — ► 0. The same 
result holds for both rigid and "free" BCs, 8u' r (r\) = 6u' r (r2) = (these conditions differ 
from ours in the respect that the configuration is bounded by a complete vacuum; i.e. 
one devoid of external fields.) Thus it would appear that one needs an infinite to 
stabilize the system. Our results are clearly at odds with DK in this respect. Although 
the finite-sized stable region found by DK was also found here, we have shown that such 
a region exists only in the presence of rigid boundaries, and then only for a < b. 
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Kumar, Coleman, <fc Kley (1994) (KCK) concluded, on the basis of the sufficient 
stability criterion (1.2), that "toroidal fields only destabilize the flow". As regards 
the VC instability, we have found that the opposite is in fact the case, at least when 
we consider the "principal range" 3/2 < a < 2, b > 1. It is only in the large-field 
(Y<t> ^ r o^o) regime that B<p destabilizes. Had the authors continued their thin-shell 
calculation to O(V^), they would have discovered that the correction to oj at peak 
growth is 

<^2,max = ^(b- a 2 /2 + a 3 /8), 

which is always damping provided that 3/2 < a < 2 and b > 1. 

As regards the enhancement of the instability for free boundaries, we note that the 
global energy change due to the perturbations, 6£, consists of three different contribu- 
tions, in general. The first is the energy change in the fluid interior, derived by KCK 
as 

6£ F = 7T J ^^p- - J"<5B x C + 2prQQ'|6-| 2 ^ rdr, (6.1) 

where £ = Su/iu + rQ'6u r (j) is the Lagrangian displacement vector. The second contri- 
bution is due to perturbations of the external vacuum field, 

6£ v = j[ \SB\ 2 rdr, (6.2) 

4 J vacuum 

while the third is a surface contribution, 6£s, which vanishes unless the equilibrium has 
surface currents (cf., Schmidt 1966). We avoid the latter here, and so the effect of free 
boundaries is given entirely by the integral (6.2), which is always positive. This led KCK 
to conclude that "stability criteria are not affected" by the BCs. However, one should 
be careful upon drawing such a conclusion from sufficient, but not necessary, criteria. 
In fact, as noted by Bateman (1978), there are numerous instances when free boundary 
instabilities grow faster than fixed boundary ones, even though 6£y > 0. The reason for 
this is simply that by allowing £ ^ at the edge of the fluid, free-boundary instabilities 
can make more effective use of the internal fluid potential energy, represented by the 
first two terms on the RHS of equation (6.1). We found ample evidence of such behavior 
in the preceding sections, and in CPS. 

Blaes & Balbus (1994) (BB) considered two-fluid models of ions and neutrals coupled 
by collisions, ionization, and recombinations. Their analysis is local, but includes an 
equilibrium azimuthal field. They found that B§ can alter the stability criterion only 
in the limit of ionization equilibrium (as opposed to ion conservation), and can in fact 
produce total stabilization for B§ J> WB Z if the ion-neutral collision frequency is below 
a certain threshold. In all other cases, B^ can cause a small reduction in growth rate, 
but does not affect the stability criterion (i.e. the critical Alfven frequency for stability 
is unchanged from the B^ = case) 6 . They take c s = 10V Z , so that the critical V<p for 



A point of formalism is worth stressing here. The finding that B,p does not affect the stability criterion, regardless 



- 29 - 



stability is Y$ ~ c s . This differs from our result, V^^rit ~ r 0^0 5 since BB's compressible 
model is sensitive to the coupling between magnetosonic and rotation-modified Alfven 
modes, which is stabilizing. BB's model does not include vertical gravity, however, so 
buoyancy instabilities which would be expected to become important near ~ c s were 
not detected. 

Gammie & Balbus (1994) (GB) considered an accretion disk model which was local 
in the radial coordinate, but global in z\ i.e. they solved for the vertical eigenmodes. 
One should be cautious in comparing our results directly to theirs, but their vertical 
node number n should compare roughly with our k, and their radial wavenumber k 
with our radial node number n. The near-coincidence of notation here is unfortunate; 
let us unambiguously re-label these parameters as n z , k z , k r , and n r , respectively. For 
a Keplerian disk, they plotted curves of constant growth rate in the (V^, V z ) plane for 
k r = and n z = 1 (their Fig. 2), finding that stabilization is achieved for V z ~ 1.5 
irrespective of V§\ this value agrees quite well with the free-boundary results of CPS 
(we found V ZjCr it — 1-43 for a = 1.5). Their BCs are similar to ours in the sense that 
far from the disk, the field lines move about freely, exerting no stress on the disk. 

On the other hand, although GB find that the growth rate decreases for increasing 
V(j, (they consider values up to V§jc s = 5), it apparently never vanishes, nor does 
affect the stability criterion. The discrepancy between these results and those of the 
present paper could be telling us something about the relative importance of vertical 
motions (which they treat in detail, and we do not) and radial ones (vice- versa). To 
date, nonlinear calculations of the BH instability have indicated that inward and out- 
ward radial motions at different z (the so-called "channel solutions" ) are the immediate 
outcome of the linear stage of the instability. It may be that the unstable modes are 
more sensitive to variations in radial structure than in vertical. GB's local approxima- 
tion in r could therefore have missed the most important effect associated with strong 
B^; namely, the prevention of the channel solution from ever forming. 

Due to the apparent similarity between GB's Fig. 2 and our Fig. 6, one might 
be tempted to make a direct comparison between the two. We caution the reader 
against it, for the following reason. The results in the former figure are for the longest 
vertical wavelength (n z = 1 or k z = 0) mode only. For this mode, the V z — ► limit 
is automatically stable, since ojj ~ — ► 0. By contrast, our critical stability curves 
are mode-independent] i.e. they reflect the requirements for stability to perturbations 
of arbitrary k. This explains the rather puzzling feature of GB's Fig. 2 in the — ► 
limit, namely, that the absolute maximum growth rate is attained not at V z = as in 
CPS, but at V z ~ 0.85. As an example, consider the a = b = 2 case, whose growth rate 
is given by the imaginary part of equation (4.5). In the k — y limit, E — ► E2(1)/k 2 , 



of its strength, is not surprising in a purely local model such as that of BB. This is because terms behaving as B^/r, 
which are crucial in the global model, are ignored in local calculations. The disappearance of B^ from the stability 
criterion in the latter case can be seen immediately from equations (4.7) and (4.17). 
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giving 

oji = -k(Vl + E 2 V 2 - 2E 1 2 /2 V z ) 1 / 2 /E 2 . 

Considered as a function of V z , the maximum of ojj occurs at V z = 0.52, independent 
of V^. The point this argument overlooks is that as V z becomes small, k necessarily 
becomes large for the most unstable mode; e.g., when V z = 0.05 and J> 0.7, there 
are no unstable modes whatsoever at k = (Fig. 3b). GB's Fig. 4 in fact shows that 
n z = 1 is not the fastest growing mode for nonzero V^. One should therefore not treat 
GB's Fig. 2 as our Fig. 6; i.e. as a critical stability diagram. 

Finally, we note that in the context of uniformly rotating magnetic stars, which 
are expected to have distributions of Ba, increasing with radius, Pitts <fc Taylor (1985) 
identified an instability having the same characteristics as the LFI (i.e. stability was 
ensured for low m (azimuthal wavenumber) modes provided that roOo J> V^o), but did 
not obtain detailed growth rates or critical stability curves. 

6.2. The Large-Field Instability: Possible Environments 

The results of §4.6. suggest that the LFI is not likely to be a threat in standard thin 
accretion disks. In some environments, however, the characteristic value for the LFI, 
V^/tqVLq J> 1, might in fact be achieved. Recent observations of flattened structures in 
massive star-forming regions (e.g., Aitken et al. 1993) suggest that such 'pseudo-discs' 
are very massive (~ 10 3 Mq) and also that the dominant magnetic field component is 
toroidal. Such massive objects are likely to be self-gravitating and sub-Keplerian, so 
that rotation may not be as important a mechanism of support as in thin disks. It 
remains to be seen how the LFI is affected by self-gravity 

On larger scales, roughly 50 % of giant molecular clouds and somewhat fewer indi- 
vidual dark clouds and cores (Goldsmith <fc Arquilla 1985) possess measured velocity 
gradients which have been interpreted as being due, at least in part, to large-scale ro- 
tation (Blitz 1993). As the magnetic fields in such objects are substantial (magnetic 
energy ~ gravitational energy ~ kinetic (nonthermal) energy; cf., Myers & Goodman 
1988), the condition V^/tqVLq J> 1 is likely to be satisfied in at least some regions. Of 
course, the effects of compressibility and self-gravity are also likely to be important, so 
a new model is needed. 

A concrete example displaying appropriate conditions for the LFI may already exist. 
The L1641 region of Orion A consists of several low-density filaments, whose major axes 
run in a roughly north-south direction. In addition to a north-south velocity gradient 
which extends across all of Orion A (~ 8 km s _1 ), L1641 also contains an east-west 
gradient, ~ 2 km s _1 , indicating that the overall velocity field of Orion A is helical 
in nature (Bally 1989). Further, the surrounding magnetic field displays the same 
symmetry (Heiles 1987). It is well-known that such a helical field is characteristic of 
superposed vertical and azimuthal fields. While figures for L1641 alone are hard to 
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come by, the average east-west gradient in the Orion A cloud as a whole (40 pc X 2 
pc) has been estimated at 0.135 km s _1 pc -1 (Kutner et al. 1977, Genzel <fc Stutzki 
1989). If this is entirely due to rotation of a cylindrical region ~ 20 pc in radius, then 
a crude estimate of the rotation velocity gives V c ~ 2.7 km s _1 . Comparing this with 
Va — 1-8 km s _1 , the density-averaged Alfven speed for the region (Heiles et al. 1993), 
one obtains Va/V c — 0.67. Given the likelihood that Va ~ (due to the predominantly 
toroidal appearance of the field), and that V c is probably a smaller contribution to the 
overall shear, one sees that values of V^,/(roOo) J> 1 should not be out of reach in this 
environment, and perhaps several others. 

6.3. Summary 

In this paper we have examined a variety of magnetic field distributions and ori- 
entations, with the principal intent of gauging their effect on the VC instability of 
magnetized accretion disks. The main results are: (1) An azimuthal field, varying as 
some inverse power of radius, has a stabilizing effect on the VC instability if its char- 
acteristic Alfven speed, V^o, is less than the characteristic rotational speed, roOo- (2) 
If Vaq J> ro^o, the system is susceptible to the LFI, whose peak growth rate increases 
with V^. This instability is more likely to affect thick, massive disks and molecular 
clouds than thin accretion disks. (3) Our calculations for finite critical wavenumbers 
suggest that complete stabilization of thin disks by an equilibrium is unlikely, since 
the required field (V<p ~ t^Qq ^> c s ) is prone to Parker instability. (4) In contrast to 
CPS, taking account of the disk's free boundaries gives qualitatively different behavior. 
In particular, whereas absolute stability can be achieved for certain rigidly-bounded 
configurations, none of the freely-bounded equilibria we examined are similarly stable. 
(5) In the absence of an equilibrium azimuthal field, a disk with a radially- varying ver- 
tical field has a smaller VC growth rate than in the constant field case. However, the 
most unstable wavenumber for fields which decrease extraordinarily quickly with radius 
may be unaffected or even increased. 

The advantages of adopting a global analysis to address questions of stability in the 
presence of strong magnetic fields are even more apparent in the present work than in 
CPS. In particular, our results show that differentially rotating gaseous bodies threaded 
by strong azimuthal, but weak vertical fields should be highly unstable for certain spe- 
cific rotational and azimuthal field profiles (§4.5.2.), a result not definitively shown by 
any local or thin shell analysis. It is hoped that future work will focus on these particu- 
lar profiles, in order to more fully examine the consequences of the ensuing instabilities. 

We thank Peter Sutherland for reading an earlier version of the manuscript and 
Omer Blaes for several useful discussions. C.C. is grateful to McMaster University for 
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APPENDICES 



A. THE PERTURBATION EQUATIONS 



We begin with the equations of ideal MHD in cylindrical polar coordinates (r, (f>, z): 



P 



dt 



+ (u-V)u 



= -pV* - v j p + B — ^) + -^(B-V)B, 

07T / 47T 



„ , 

— = Vx(uxB). 



V-u 



V B = 0. 



(A.l) 

(A.2) 
(A.3) 



Here p is the gas pressure, p the constant density, u the fluid velocity, and \I/ = —GM/r 
the gravitational potential. Substituting perturbations of the form (3.1) into these 
equations and only retaining terms of linear order in perturbed quantities, we obtain 



ikB ( Ba 

iujSvl + V6h -SB + 2 ( -^SBa, - QSi 



4ir pr 



B \ - B' 

+2 6B r - B6u r (f> z -6B r z = 0, 

\ 4irp I 4irp 

iuSB — ikB z 6u + 2(ASB r - A6u r )(j> + B' z 6u r i = 0, 
— (rSu r )' + ikSu z = 0, 



(A.4) 
(A.5) 
(A.6) 



where h = p/p + B 2 /(8irp) is the specific enthalpy, A = -rO'/2, B = -[(rO)' + 0]/2] 
are the usual Oort shear parameters, and A = — r(_B^/r)'/2, B = — (B 1 ^ + B^jr)j2 
are their magnetic counterparts. For future reference, we note that with the power-law 
forms (2.1), the first two of these equations become 



ikB / Ba 

iujbxi + VSh -SB + 2 ( -^SBa, - Q6u, 



+ 



4-Kp 

(2 - a)nSu r - ^~ b)B4, 6B.* 



Ait pr 



iujSB — ikB z 6u + I aQSB 



bB, 



Ait pr 
(1 - c)B . 
Aitpr 

' + (1 — c) — -6u r z = 0. 



z 6B r z = 0, 



(A.7) 
(A.8) 



Note that equations (A.5) and (A.6) imply V ■ <5B = 0. Resolving equations (A.7) and 
(A.8) into components, and using equation (A.6), one can eliminate all variables except 
Su ri leading to the perturbation equation (3.2). 
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B. ROTATING VS. NONROTATING EQUILIBRIA 



The stationary pressure distribution in the nonrotating case can be found from 
equation (2.2) with Oo = 0. When V z = constant, the pressure maximum relation is 
simply 



GM 
ro 



(B.l) 



Note that this requires 6 > 2 for a sensible equilibrium. Using this to eliminate GM in 
equation (2.2) and integrating gives 



n=o 



ro 1 - (r/rp)- 2 ^ 1 ) 

r 2(6-1) 



(B.2) 



For a constant vertical field, equation (2.5) reads (in proper units) 



p 
P 



— + (n,n ) 2 x 




(r/rp)- 2 ^- 1 ) 
2(a - 1) 



(2-6) 



VI 



r 2 2 

1 o li o 



r 



(r/r ) 



-2(6-1) 



When 6 = a, this becomes 
P Po 



P 



P 
Po 



+ [(r O ) 2 - (2 - b)Vjo] 



1 + 



1 - (r/r ) 



2(6-1) 



-2(6-1) 



= j + Vj (b eJf -2) 



^-1 + 

r 



where 



b eff = 



2(6-1) 
1 _ ( r / ro )-2(fc-l)' 



+ 6. 



(B.3) 



Clearly, equation (B.3) is identical to equation (B.2), but for the replacement of 6 
by beff - The two equations become identical in the limit V<po ^> roO(), which is precisely 
the regime of the LFI found in this paper. In fact, as soon as Vqq J> roOo, one would 
expect that the rotating system should start to display much of the qualitative behavior 
of its nonrotating counterpart, since then the contribution of the magnetic terms to the 
pressure is of the same sign in equations (B.2) and (B.3) for 6 J> 1. Finally, one might 
be tempted to blame the LFI for 6 < 2 entirely on the violated equilibrium condition 
(B.l). However, this condition applies only when V^o ^> tqQ,q. As an example, take 
6 = 1.7. Then Fig. lb shows that all equilibria with < Va, < 1.83 are allowed. 



C. THIN-SHELL APPROXIMATION 
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Following Kumar, Coleman, <fc Kley (1994), we adopt a thin shell approximation, 
in which the radial dependence of equilibrium quantities is ignored to first order, but 
their derivatives are not. The perturbation equation (4.1) then becomes 

^ + « = o, (C.l) 

where 

Q = ^(bu 2 - 20i) + ~ auj 2 ) - _ {k 2 + 3/4)j (C . 2) 

UJ UJ UJ 

i\) = y/r6u r , and ( = r — 1. Since Q$ is a constant, the solution of equation (C.l) is 
%jj = c\ sin y/Qo C + c 2 cos y/Qo Ci applying the rigid BCs then gives Qq = (nir/s) 2 , 
where s is the shell half-thickness and n the radial mode number. Assuming ks n, 
equation (C.2) yields the following characteristic polynomial: 

+ [2(a - 2 - bV}) - 2Q 2 a ]uj 2 + 8n A V<puj + 0^0^ + 2(6 - 2)V$ - 2a] = 0, (C.3) 

where the reader is reminded that all equilibrium quantities are to be evaluated at 
r = tq = 1. 

The roots of equation (C.3), although calculable analytically, are algebraically com- 
plicated and do not give much physical insight. Kumar et al. (1994) adopted a procedure 
equivalent to expanding uj in powers of V^, taking the latter as a small quantity. As our 
object is to obtain the critical stability curves, it is more useful for our needs to place 
no restriction on V§\ rather, we take as a small parameter. Expanding uj as 

UJ = UJ + O y 46D 1 + + ■ ■ ■ , 

substituting into equation (C.3), and solving the resulting equation in orders of Q^, we 
find there are two branches of the dispersion relation. One gives all real contributions 
to uj; the other has ujq = and 



-1V$ ± [b(b - 2)TA 4 + 2(a + b - ab)V 2 - a{2 



11/2 



a - 2 - bV 2 

Positivity of the square-root argument leads to the stability criterion (4.14). 

D. PROOF THAT uj 2 IS REAL WHEN B z = B z (r) AND B^ = 
The perturbation equation in this case is [equation (3.2)]: 



(C.4) 



^[ruj 2 (Su r )']' + lk 2 r 



(0 2y _ o£y 



+ 4k -uj 2 (k 2 + \)} fiv r = 0. 

uj z V r z 
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Multiplying through by rSu* (an asterisk denotes the complex conjugate) and integrat- 
ing, one finds 



rr 



rri 
>r\ 

where 



k r 



(0 



2v ( v ?y 



+ 



4k 2 oj 2 n 2 



C'J 



6j' z [k 2 + \] }r\6u r \' z dr + I = 0, (D.l) 



rri r „ rri 

1= 6u;[rCo 2 (6u r )']'dr = rCo 2 (6u r )'6u* r - rCo 2 \(6u r )'\ 2 dr. (D.2) 

The latter result is obtained via integration by parts. 

We consider rigid BCs only, since we restrict consideration to B z = constant in 
the free boundary case to avoid currents (§2.2.1.). Applying 6u r (r\) = 6u r (r2) = to 
equation (D.2), the first RHS term vanishes. Substituting the result back into equation 
(D.l) and taking the imaginary part of the entire expression gives 



(0J 2 )j 



r2 
n 



\(6u, 



\/|2 



+ 



C'J 



212 z 



|^u r | 2 > rdr = 0, 



where a subscript I indicates the imaginary part. The integrand is positive definite for 
all r, showing that (oo 2 )i = 0. 
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